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1. Introduction 

Fractional-Order systems, or systems containing fractional derivatives and integrals, have been studied by 
many in the engineering area (Heaviside, 1922; Bush, 1929; Goldman, 1949; Holbrook, 1966; Starkey, 
1954; Carslaw and Jeager, 1948; Scott, 1955; and Mikusinski, 1959). Additionally, very readable 
discussions, devoted specifically to the subject, are presented by Oldham and Spanier (1974), Miller and 
Ross (1993), and Podlubny (1999a). It should be noted that there are a growing number of physical 
systems whose behavior can be compactly described using fractional system theory. Of specific interest to 
electrical engineers are long lines (Heaviside, 1922), electrochemical processes (Ichise, Nagayanagi, and 
Kojima, 1971; Sun, Onaral, and Tsao, 1984), dielectric polarization (Sun, Abdelwahab, and Onaral, 
1984), colored noise (Mandelbrot, 1967), viscoelastic materials (Bagley and Calico, 1991; Koeller, 1984; 
Koeller, 1986; Skaar, Michel, and Miller, 1988), and chaos (Hartley, Lorenzo, and Qammar, 1995). 

With the growing number of applications, it is important to establish a theory of control for these 
fractional-order systems, and for the potential use of fractional-order systems as feedback compensators. 
This topic is addressed in this paper. The first section discusses the control of fractional-order systems 
using a vector space representation, where initialization is included in the discussion. It should be noted 
that Bagley and Calico (1991) and Padovan and Sawicki (1997) both present a fractional state-space 
representation, which do not include the important historic effects. Incorporation of these effects based on 
the initialized fractional calculus (Lorenzo and Hartley (2000)) is presented here. The control methods 
presented in this paper are based on the initialized fractional order system theory (Hartley and Lorenzo 
(1998)). The second section presents an input-output approach. Some of the problems encountered in 
these sections are a) the need to introduce a new complex plane to study the dynamics of fractional-order 
systems, b) the need to properly define the Laplace transform of the fractional derivative, and c) the 
proper inclusion of the initialization response in the system and control formulation. Following this, the 
next section generalizes the proportional-plus-integral-control (Pl-control) and PID-control (PI-plus- 
derivative) concepts using fractional integrals. This is then further generalized using general fractional- 
order compensators. Finally the compensator concept is generalized by the use of a continuum of 
fractions in the compensator via the concept of order-distributions. The last section introduces fractional 
feedback in discrete-time. 
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Vector Space Methods for Fractional-Order Control 


2.1 General Vector Representation and Initialization 

A useful representation for systems of fractional-order differential equations is the vector space of 
fractional dynamic variables. This representation is a generalization of the state equations of usual 
integer-order system theory. It will be important to notice that for fractional-order systems, the “state” of 
the system is not given by the dynamic variable vector, because the initialization vector, traditionally a 
vector of constants, has been shown to be time varying [Lorenzo and Hartley (2000)]. This section 
summarizes the multivariable initialized fractional -order system theory. 

For a chosen q, typically 0<^<1 which implies a fractional-order system, the vector 
representation can be written as 

.. Dj x(t) = /Ly(t) + Bu(r) , y/(x,q,a,c,t ) given for t>c (2.1.1) 

y(r) = Cx(t) + Du(t) . (2.1.2) 

where x(t) , y(t) , u(t) , and \ff(x,q,a,0,t) can generally be considered to be vectors of the appropriate 

dimension, with A , B , C, and D now being matrices of the appropriate dimension, and the starting time, 
c, will usually be chosen as zero. The time period a<t<c is called the initialization period. It is 
assumed that .v(r)=0 for all t<a . The choice of q is problem dependent and is discussed further below. 
The derivative notation is 


0 O/ x(t) = 0 d x(t) + \jf(x,q,a,0,t) 


(2.1.3) 


where () d? is the un-initialized fractional derivative starting at r = 0, and , the vector 
initialization function, is determined as 


\j/(x,q,a,Q,t) = a d q a x{t) = — 
~ dt 

or more specifically we mean 


1 f X(T) 

n\-q)[( t - T y 


dz 


W-, (x,,q,a,0,t) s a d q x,(t) = — 

at 


Xj (r) 

ni-«)](r-j -r 


1 u 

—I 

-nW 


dr 


(2.1.4a) 

(2.1.4b) 


Note that only initialization through terminal initialization is considered in this paper, although more 
general initializations are discussed in Lorenzo and Hartley ( 1 998). 

Equation (2.1.1) can be rewritten using Equation (2.1.3) as 


0 d q x(T) + yf(x(t),q,a,0,t) = Ax(t ) -I- Bu(t) , 


(2.1.5) 
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where the initialization function is the vector of functions. 


i//(x, q, a,0,r ) = 


y/(x { ,q<a,0, t ) 
yi(x 2 ,q,a,0,t) 

y/(x n ,q,a,0,t) 


In what follows, this will, without loss of generality, be written as 


V(x,q,a,0,t)= V >_(t) 


¥\ (0 
y/ 2 (t) 


¥n(t) 


( 2 . 1 . 6 ) 


(2.1.7) 


An alternative representation that can be useful for determining the initialization response is to 
rewrite the system dynamic equations into the fractional integral form, 

x(t) = 0 D; q [Ax{t) + Bu(t)] (2.1.8) 

which can then be written as 

x(r) = 0 d~ q [Ax(r) + B«(r)]+ yr(Ax(t) + Bu(t),-q,a,0,t ) . (2.1.9) 

This formulation allows a somewhat simplified determination of the initialization response 
V f(Ax(t) + Bu(t),-q,a,0,t) with respect to Equation (2.1. 1), as only a fractional integral of the history of 
x(t) is required here. The history of x(t) is easily obtained for a given input in negative time. 

The fractional dynamic variables in the system of Equations (2.1.1) and (2.1.2) are not states in 
the true sense of the name “state” space. In the usual integer-order system theory, the set of states of the 
system, known at any given point in time, along with the system equations, are sufficient to predict the 
response of the system both forward or backward in time. Otherwise stated, for integer-order systems, the 
collection of numbers x(t ) , at any time t , specify the complete “state” of the system at that time, and 
with the differential equations, for all future times. Mathematically stated, for integer-order systems, the 
system will have a unique time response given its initial state. It should be clear, however, that the 
fractional dynamic variables do not represent the “state” of a system at any given time due to the 
requirement of the initialization function vector, which carries information about the history of the 
elements of the system. Consequently, as the initialization function vector is generally required, the set of 
elements of the vector x{t ), evaluated at any point in time, does not specify the entire “state” of the 
system. Thus for fractional-order systems, the ability to predict the future response of a system requires 
the set of fractional differential equations along with their initialization functions, that is. Equation 
( 2 . 1 . 1 ). 

The availability of the fractional vector formulation allows many choices for the basis q . Thus, the 
total number of fractional dynamic variables can increase if the basis q is chosen to be smaller. The least 
number of fractional dynamic variables is obtained by choosing the basis q as the largest common 
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fraction of the differential orders. For example, a system with two elements of order q l =1/2, and 
<72 =1/3 , would require a basis q of at most q = 1/6 . Clearly, such a system could also use a basis q of 
1/12, 1/18, 1/24, or any other fraction of the form 1/6/?, «= 1,2,3,. . . . If the original differential equations 
have non-rationally related derivatives, such as 0 d' ly ^x(t) 4- n d) l!t x(t) + x(t) = u(t), then this 
approach will require an approximation of the value of q , as integer sub-multiples of the orders are not 
obtainable. 

As in the case with integer-order systems, a particular input-output representation can have an 
infinite number of vector space representations. We can change vector variables in Equations (2.1.1) and 
( 2 . 1 . 2 ) using 


x(t) = Tz(t) (2.1.10) 

where T is a square matrix of the same size as the vectors x(t) and z(t ) . With this definition, Equations 


( 2 . 1 . 1 ) and ( 2 . 1 . 2 ) become 

0 d? Tzit) + V(Tz(t\q,a$j) = AT z(t) + Bu(t ) (2.1.12) 

y(t) = CTz(t) + Du{t). (2.1.13) 

Rearranging the first of these gives the set to be 

o d?z(t) + T- l y/(Tz(t),q,a,0,t ) = T~ l ATz{t) + T~ l Bu(t) (2.1.14) 

y(t) = CTzit) + Du(t) (2.1.15) 


Thus we have changed the set of vector variables from ^(/) to z(.t). Notice that the new system 

matrix is T~ l AT , the new input matrix is T~ l B , the new output matrix is CT , and the direct feedthrough 
matrix is unchanged as D . The w-plane eigenvalues of the system matrix are the system poles, and these 
must be unchanged with the change of variable. Thus the original system matrix A and the new system 

matrix T~ l AT must have the same w-plane eigenvalues. Thus, as in integer order systems, the system 
eigenvalues are unchanged via the similarity transformation. 

The initialization problem can be solved for A'(/) by using Laplace transforms generalized by the 
proper inclusion of the vector initialization function. As is shown in Lorenzo and Hartley (1998), the 
Laplace transforms of Equations (2.1.5) and (2.1.2) are 

s q Xfs) + y(s) = A2((s) + BU(s ) (2.1.16) 

Y(s) = C X(s) + DU(s) . (2.1.17) 

Rearranging Equation (2.1.16) gives 

(/^-Afes) = BU(s) -\j/(s). (2.1.18) 
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or equivalently 


X(s) = (ls q -A) ‘ BU(s) - (ls q - aJ'Vcs) . 


(2.1.19) 


Inserting Equation (2.1.19) into Equation (2.1.17) gives the Laplace transform of the system output 
response 


F(s) = {c(/s 9 -A)"' B + d}[/(5) - c(ls q -a)“V(5) • 


( 2 . 1 . 20 ) 


Using the system's impulse response, the matrix F-function (which is the generalized matrix exponential, 
see Hartley & Lorenzo (1998)), 


r, h '1 * '*■' I rfCn T = V ' {( ,s ’ - ' }• 1 > 0 • 

nn( ,i + 1 >4) 

Equation (2. 1.20) can be inverse Laplace transformed to give the overall time response 

r t 

y(t) = Cj F q [A,r]Bu(t -T)dT + Du(t) - C j F q [A,t]y(t -T)dr . 

0 0 

This can also be represented as 

y(t) = c|F ? [A,r]^w(r-r) - i//(i-t)}c/t + Du(t). 

o 

or equivalently as 

y(t) = C j F q [A,t - t]{b u(t) -y{T)}dT + Du(t). 


( 2 . 1 . 21 ) 


( 2 . 1 . 22 ) 


(2.1.23) 


(2.1.24) 


2.2 Stability Analysis in the w-plane 

To understand the dynamic behavior and stability properties of the system of Equation (2.1.1), it 
is necessary to analyze the eigenvalues of the system A-matrix. For integer-order linear system theory ( q 
= 1), the eigenvalues of the A-matrix are studied in the complex Laplace 5-plane. The stability boundary 
in the 5-plane is the imaginary axis; any poles lying to the left of the imaginary axis represent a stable 
time response, while the poles lying to the right of the imaginary axis represent an unstable time response. 
Examining Equation (2.1.19), however, indicates that the eigenvalues of the A-matrix must now be 
evaluated in what would appear to be the s 9 -plane. Rather than dealing with the fractional power of 5 , the 
analysis is simplified if a change of variables is used. We will define w=s q , and then the eigenvalue 
analysis will be performed in the new complex vv-plane, which is a mapping of the 5-plane. 
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To accomplish this, it is necessary to map the 5-plane, along with the time-domain function 
properties associated with each point, into the complex w-plane. To simplify discussion we will limit the 
order of the fractional operator to 0 < q < 1 . Let 

vv = pe J ‘ t> = a + jp and s = re i6 . (2.2.1) 

Then referring to the definition of w, 

iv = 5 ? = (re jd y = r q e jqB = pe J * (2.2.2) 

which gives 


p-r q and <p = q6 . (2.2.3) 

With this equation, it is possible to map either lines of constant radius, or lines of constant angle from the 
5-plane into the tv-plane. Of particular interest, is the image of the line of 5-plane stability (the imaginary 

axis), that is s = re ±iK '- . The image of this line in the w-plane is 

w = r*e ±J *%, (2.2.4) 

+ qK 

which is the pair of lines at 0 = — - — . Thus, the right half of the 5-plane maps into a wedge in the 
tv-plane of angle less than ± 90 q degrees, that is, the right half 5-plane maps into 


0 


qn 


(2.2.5) 


For example, with q = , the right half of the 5-plane maps into the wedge bounded by |0| < ^ . 

A half-order system with its tv-plane poles in this wedge, that is with \<f>\< 7 y / ^ , would be unstable, and the 
corresponding F-function response would grow without bound. 


It is also important to consider the mapping of the negative real 5-plane axis, s = re ±Jn . The 


image is 


vv = r q e ±iqn . (2.2.6) 

Thus the entire primary sheet of the 5-plane maps into a vv-plane wedge of angle less than 
±180<3 i degrees, while all of the secondary 5-domain sheets map into the remainder of the w-plane. 
For example, if q = y 2 , then the negative real 5-plane axis maps into the w-plane lines at ± 90 degrees. 

To summarize the above, the shape of the F-function time response, F [a,r], depends upon both 
q , and the parameter a , which is a particular eigenvalue of the system A-matrix of Equation (2.1.1). For 
a fixed value of q , the angle 0 of the parameter a, as measured from the positive real w-axis, 
determines the type of response to expect. For small angles, \<p\< qn /_ , the time response will be unstable 
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and oscillatory, corresponding to poles in the right half 5-plane (see plots in Hartley and Lorenzo (1998)). 
For larger angles, qK / 4<\(p\<q7T , the time response will be stable and oscillatory (underdamped), 

corresponding to poles in the left half 5-plane. For ]<(>) = £/zr , the time response will be overdamped. For 
even larger angles, |0| > <77T , the time response will be hyperdamped, corresponding to poles on 
secondary 7 Riemann sheets. 

All of the usual control techniques concerning eigenvalues, or poles, can be used in the w-plane. 
As is discussed in what follows, pole placement using state feedback will place the closed-loop poles in 
the w-plane, and all the root-locus analysis is performed in the w-plane as well. 


2.3 Fractional Vector Feedback 

This section considers the use of fractional vector feedback for control design. The system is the 
vector space of fractional -order elements from Equations (2.1.1) and (2. 1.2) and has the form 

0 d? x(t)+ y/(x,q,a,0,t) = Ax(t) + Bu(t), y(t) = C x(t) + Du(t) . (2.3.1) 

Typically, vector feedback is implemented in the form 

u(t) = - K x(t) + r(t) (2.3.2) 

where r is a vector reference input, and K is the feedback gain matrix to be determined. The closed-loop 
system then becomes 

0 d? x(t) - [A-BK]x(t) - \y(x{t),q,a,0,t) + Br(t) , y(t) = [C - DK]x(t) + D r(t) (2.3.3) 

By choosing K appropriately and using any standard pole placement tools, such as the Bass-Gura 
method or Ackerman's method (Kailath (1980)), it is possible to place the system eigenvalues anywhere 
in the w-plane. It should be remembered that the Bass-Gura approach is based on similarity 
transformations of the vector space. As shown above, similarity transformations can be applied directly to 
fractional-order systems given in vector space, exactly as in the integer-order situation. 

A particular problem associated with fractional-order systems is the presence of the initialization 
term on the right side of Equation (2.3.3). In the usual integer-order equation (<7=1), the initialization term 
is the familiar initial condition vector for all the states. Now, however, the initialization function is no 
longer a constant function of time, but is a time varying function. It is important to realize that the 
initialization term will always be a decaying function of time (for 0«?<1). To understand this, it will be 
assumed that the initialization begins at some time a less than zero but not equal to negative infinity, that 
the input function u(t) is of finite size except for the possibility of delta functions, and that the order of the 
fractional derivative is 0<g< 1 . With these assumptions, and remembering that the initialization function is 
the future response of only the fractional derivative term due to all past inputs to the system, it must 
necessarily decay to zero. This can be seen from the fact that for 0<q< 1 , the fractional derivative is 
effectively a lossy element, and even if the system response is unstable, the finite negative a assumption 
requires that the energy put into the fractional element is finite. Thus the energy coming back from the 
fractional element in positive time is necessarily finite. Thus we may proceed with the implementation of 
vector feedback controllers without the worry of an unstable plant in negative time generating an unstable 
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initialization response in positive time, implying that any vector feedback controller would never be able 
to stabilize the response. 

It should also be noted that a standard vector feedback approach, the LQR design method, does 
not have any particular meaning for the fractional-order vector space system at this time. It can, however, 
still be used directly, using A and B and any positive semi-definite Q and R matrices, to give a guaranteed 
stable, hyperdamped system. 


2.4 Observers for Fractional-order Systems 

Just as it is in integer-order system theory, it is important to create observers, or vector estimators, 
for fractional-order systems. This section will present the theory necessary for designing fractional-order 
system observers. The fractional-order vector estimator has the form 

0 d? x(t) + \fi(x(t),q,a,0,t) = Ax{t) + Bu(t) - L^(r)-y(r)], y(t) = C x(t) + D u . (2.4.1) 

where a non-zero initialization function, rfr , has been assumed for the observer. The vector error e(t) is 
defined as the difference between the real system output *(?) , and the estimated observer output x(t) 


e{t) = x(t)-x(t) . 


(2.4.2) 


The observer error gain L is determined so as to force the error between the two plant vectors to go to 
zero. The dynamics of the error are obtained by fractionally differentiating Equation (2.4.2), 

0 d?e(t)= a d?x(t)- 0 d?x(t) (2.4.3) 

Substituting the system equations from Equations (2.4.1) and (2.3.1) yields 

a df e(t) = [4x(r)+ Bu(t)-\fHx(t),q,aS)J)\- [a*(0 + Bu{t)~ \jf(x(t),q,a,0,t) - L[y(t)~ y(f)]]. (2.4.4) 

Now replacing the sensed system outputs, y(t) and y(r ) , with the vector variables using Equations 
(2.3.1) and (2.4.1), yields 

0 d?e(t) = [A4r)+fiw(r)-i//(^,^,a,0,r)j- 

r — ~ -i (2.4.5) 

\Ax(t)+ Bu(t)~ ty r (x(t),q,a,Qj) - L[(C ; t(r) + Dw(r))-(C ; v(r) + D«(r))]J 

Eliminating the terms that subtract out, replacing e{t) = x(t)~x(t ) , and combing terms, gives 

0 d?e(t) = (. A-LC)e{t ) - ( y(x(t),q,a,0,t ) - \pr(x(t),q,a,0,t)). (2.4.6) 

The matrix L is determined to force the observer error to zero by placing the eigenvalues of A-LC in a 
stable region of the w-plane using standard methods. As discussed previously, the initialization response 
eventually decays to zero for any system for 0<q< 1 , and only has a transient effect on the observation 
error, however a proper choice of i)/ will allow the error to go to zero sooner than if iff was simply zero. 
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2.5 Fractional Vector Feedback with Integer-order Systems 

An advantage of fractional-order controllers for integer-order systems lies in the use of fractional 
derivatives for feedback purposes. As an example, a simple scalar plant is given by 

i(/) = -x{t) + u(t) , y(r)=x(r). (2.5.1) 

This can be described, however, using fractional-order vector elements. Choosing < 7 = 1/2, this system 
becomes 


0 d° 5 x,(r) = x 2 (t) 

0 d? 5 x 2 (t) = -V, it) + u(t) (2.5.2) 


y(t) = .v l (f) . 

Fractional vector feedback can now be used, as discussed in the previous section, to place the poles 
wherever desired in the u-plane for < 7 = 1/2. For a vector feedback gain vector K, the closed loop 
characteristic equation is 


A(w) = \wI-A + BK\ = w~ + k 2 w+(l + k l ) = 0- (2.5.3) 

Then k 2 and k { are chosen to allow desired placement of the closed-loop poles. Assuming that we can 
measure the “substate” x 2 , we gain the ability to obtain underdamped responses in this first-order system 
by using negative values of the control gains. Even more freedom of response design can be obtained by 
measuring any arbitrary number of fractional substates. Thus this simple system could be written, with 
q=l/N, N = 2, 3, 4, ... as 


0 d? x x (t) = x 2 (t) 

0 d? x 2 (t) = x 3 (t) 

: (2-5.4) 

0 d? x N _\(t) = x N _ 2 (t) 

0 d? x N (t) = x,(r) + u(t) 

y{t) = x { (t) 

Now by using fractional vector feedback, the closed-loop characteristic equation in the w-plane with 
q=l/N , is 

A(w) = \wI-A + BK\ = 0 , 
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or 

A(vv’) = vi’ + A' jV v v‘ v 1 +A’ a ,_jH’" v ” + +A'2W+(1 + A'i ) = 0 . (2.5.6) 

This allows general system response design by placement of the closed-loop poles in this vv-plane. If the 
fractional variables cannot be measured directly, then we can revert to using an observer to generate them 
as discussed in the previous section. 


3. Input-Output Methods for Control 

3.1 Essence of the Approach 

This section presents a discussion of fractional-order systems from a classical input-output 
control approach. The important results here are that the classical control tools, frequency response and 
root locus, are still applicable with the appropriate modifications. The frequency response approach 
applies directly to fractional-order systems as long as the primary roots are used in evaluating the 
individual fractional elements. Likewise, the root locus approach applies directly to fractional-order 
systems as long as the root locus analysis is performed in the w-plane. Throughout the discussion, the 
time-varying initialization function is considered to be a disturbance entering the system, and is 
accommodated with disturbance rejection techniques. 


3.2 Conversion of Vector Equations to Transfer Matrix Form 

Input-output representations can be obtained in several ways, and they provide significant insight 
into system behavior that compliment the vector space approach. The term input-output representation 
usually implies a transfer function, or transfer matrix, but we will include the initialization response for 
completeness. One way to obtain an input-output representation is to replace any physical system 
elements by their Laplace domain impedance equivalents, and then do the necessary algebra to obtain a 
transfer function. This approach works equally well with both integer-order systems and systems 
containing fractional-order elements. 

The other approach to obtaining an input-output representation is to convert an already existing 
vector space representation into input-output form. Fortunately, this is easily done in a manner similar to 
that in integer-order system theory. As presented in Equation (2. 1 .20), the conversion is 

Y(s) = jc(/5 ? -Af' B + d\u_(s ) - C (ls q -A)'V(s), (3.2.1) 

where the system transfer matrix is 

G(5) = ^’(/s ? -a)' 1 S + d} (3.2.2) 

and 

Y{s) = G(s)U(s) - c(ls q - aYv(s) . (3.2.3) 

It is just as easy to go from a transfer function representation containing fractional elements to a vector 
space representation using any of the standard canonical forms from integer-order theory, and a chosen 
value of q (see Kailath ( 1980) for a discussion of canonical forms). 
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3.3 Sinusoidal Response of Fractional-order Operators in the Time Domain 

This section presents the time response of a fractional order system to a sinusoidal input. The 
results will be used in the next section to clarify the multivalued nature of the frequency response of 
fractional operators. It should be remembered that the frequency domain approach assumes that the time 
responses are in sinusoidal steady state. Whenever an input is applied to a system, the response will 
always consist of a transient part plus a steady state part. The frequency response approach assumes that 
the transient has decayed away, and that the response is in sinusoidal steady state. For fractional-order 
systems, sinusoidal steady state also implies that the initialization response has decayed to near zero. 

Oldham and Spanier (pp. 108-1 10) give the fractional differintegral of a periodic function, and we will 
follow that discussion. We can write any periodic function with period T, as follows, 


V 1 („ „j2nktn . „ -j2nkt/T\ 

/(?) = Zj\ c k e + Cke 1 ), 

k = 1 

where the coefficients c k can be determined from the corresponding Fourier integrals 

q. =±]nt)e- J **' T dt 
0 


(3.3.1) 


(3.3.2) 


and Ck is the complex conjugate of c k . The fractional differintegral of f(t) then requires the 
differintegral of each of the separate exponentials. Oldham and Spanier show that 



(3.3.3) 


where y * is the incomplete gamma function. These two differintegrals contain both the transient and 
steady state responses. To obtain the sinusoidal steady state fractional derivative of the periodic function 
in Equation (3.3.1 ), an asymptotic expansion for large values of t of the y * terms gives 


o d? fit) 




+ Cke 



!l+i 

T 4 




V 


J 


(3.3.4) 


2 K 

Defining the radian frequency to be co 0 =— 


gives the equivalent response 


0 d^ = 


7 k(t> 0 t+~- _ -J 

c k e 1 ~ i + Cke 


k(Of)t | 1 


*=1 


(3.3.5) 
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For any given input frequency kco 0 , it can clearly be seen that the magnitude of the corresponding 
fractionally differintegrated steady-state output sinusoid has its magnitude scaled by (ka> 0 ) q , and is phase 

shifted by —■ . For example, after the decay of the transient, 0 d q sin(<yf) = (co) q sin(twr + {rcq/ 2)). This 
result generalizes the response obtained for integer-order systems. 


3.4 Sinusoidal Response of Fractional-order Operators in the Frequency Domain 

The transfer functions for integer-order systems are usually given in the Laplace domain. To 
obtain the frequency response from such transfer functions, the terms containing the Laplace variable s, 

are replaced by the radian frequency co , that is 5" —>(jco)", where n is any integer exponent. When this 
is done for fractional-order transfer functions however, the question of multiple solutions exists, as the 
substitution becomes s'* 1 -^(jco )'* 1 , where q is fractional and nq is not necessarily an integer. The 

question here is actually, which of the roots of j nq do you use, as there are generally many roots, which 
lie on the unit circle in the complex plane. The primary' root is considered to be the one with the smallest 
angle from the positive real axis, with the remaining roots being the secondary roots. The answer to this 
question is given by the time domain result of the previous section. That is, using the primary roots given 

by the frequency domain substitution s nq — > ( j(o) nq will give the frequency response corresponding to the 

correct time domain response. For example, for s 05 , substituting s=jco gives ( jm ) 05 = j 05 (O 05 . 

Recognizing that j 05 has two roots, e jn/4 and e i5n,A , the primary roots is always chosen for the 

frequency response, that is e jlz/4 . This observation then allows the use of the standard frequency domain 
analysis tools such as the Bode plot and the Nyquist plot. 


3.5 Frequency Response of Fractional-order Differintegrals 

To demonstrate the frequency domain approach, this section considers the frequency response of 
simple differintegrals. The uninitialized Laplace transform of a fractional integral (<y<0) or derivative 
( < 7 > 0 ) operation is 


L{ 0 d? f(t)}= s q f(s). 

Thus the transfer function of the operator is 

H(s) s'/W 

f(s) f(s) 


= s q . 


for q real. 


(3.5.1) 


(3.5.2) 


To obtain the frequency response, let s q —> (jco) q , thus 


H(jco) = {jco ) q . 


(3.5.3) 
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The magnitude response is simply 


\H(jco)\ = co q , (3.5.4) 

which rolls off at 20 q dB/decade on the Bode plot, and the phase shift is given by the angle of the primary 
root of j q 

ZH( jto) = ^ . (3.5.5) 

These responses are shown in Bode form in Figure 3.5.1 and in Nyquist form in Figure 3.5.2'. 



Frequency, radians 



Frequency, radians 

Figure 3.5.1. Bode frequency response for fractional-order operators, 
with q varying from -1 to + 1 in steps of 0.25. 
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Rea! H(jw) 


Figure 3.5.2. Nyquist frequency response for fractional-order operators. 

The arrow shows the direction of increasing frequency. 

3.6 Frequency Response of Ultradamped Systems 

The time domain behavior of fractional-order systems was discussed in Section 2.2 with respect 
to their u-plane pole locations. There it was shown that any poles lying to the left of the wedge, given by 
lines with angles qn in the w-plane, are on the secondary Riemann sheet of the 5-plane, and these poles 
were termed hyperdamped, as they were damped more than the usual integer-order overdamped poles. 
Now, with respect to the w-plane, it is with some necessity that we distinguish between poles that are on 
the negative real w-plane axis, which we will call ultradamped (rather than overhyperdamped), and 
complex conjugate poles that are still in the hyperdamped region, which will be called simply 
hyperdamped. 

An ultradamped system will consist of parallel combinations of systems of the form 

ytvt k 

H(s) = — - = — , u > 0 , (3.6.1) 

U(s) s q + a 

where Y is the system output, U is the system input, k is the system gain, -a is the ultradamped system 
pole in the vv-plane, and q> 0. Without loss of insight, we can let a=l, and k=\. Then we can analyze the 
frequency response of this transfer function, 
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H(j(o) = 


1 

U°>Y + 1 ' 


(3.6.2) 


For small values of the frequency, the transfer function is H( jco) = 1 , and the magnitude is thus 1, and the 

phase lag is 0. For large values of the frequency, the transfer function becomes H( j( 0 ) = (j(o)~ q , and the 
frequency response reverts to that of the simple fractional-order operator discussed in the previous 
section. Bode and Nyquist plots of this behavior are presented in Figures 3.6.1 and 3.6.2 respectively. 
With poles and zeros in these locations, it is simple to construct both phase lead and phase lag systems. 
These will be discussed later with regard to feedback to control. 


3.7 Frequency Response of Hyperdamped Systems 

Hyperdamped systems have a pair of complex conjugate poles off the negatiye real axis of the 
w-plane. Here there are several types of system behavior depending upon the specific locations of these 
poles in the w-plane. These behaviors are addressed in this section. Several definitions are given first 
along with corresponding pole locations in the w-plane. 

We will first address systems that can be realized with passive energy storage elements. A passive 
energy storage element is one that cannot return more energy to a system than was placed into the element 
by the system in the past. An active element is one that can return more energy to a system than the 
system placed into it in the past, and obviously can behave as a passive element if designed so. Typically, 
an active element will have associated with it either a large gain or a negative gain. Necessary, but not 
sufficient, conditions for a fractional -order system to be passive are that its minimal transfer function 
denominator have all positive coefficients and that all of its poles lie to the left of the w-plane stability 
boundary. Another concept traditionally associated with passivity is the positive-real concept. To be a 
positive-real system, the frequency response of its transfer function must always lie in the right-half 
Nyquist plane. Stated otherwise, the maximum phase shift of a positive -real system is bounded by plus, or 
minus, ninety degrees. Passive fractional-order systems are not necessarily positive-real (Figure 3.6.2 
shows passive system frequency responses some of which are not positive-real). 

A minimum phase system has the smallest possible phase shift for a given magnitude response. 
An implication of this for integer-order systems is that all of the system’s poles and zeros must lie in the 
left half r-plane. For fractional-order systems, the implication of being minimum phase is that all of the 
poles and zeros of a given system must lie to the left of the instability wedge in the w-plane. Considering 
all-pole systems, passive and positive-real systems are clearly minimum phase, however, a fractional- 
order system can be minimum phase without being passive or positive-real. 

A simple example of this is the transfer function 


H(s) = 


1 

5 + as 05 + 1 


(3.7.1) 


Referring to Figures 3.7.1 and 3.7.2, the w-plane is shown for Equation (3.7.1) with identification of the 
regions of passivity, positive realness, minimum phase, various damping, and stability for this particular 
system, along with corresponding step responses for various values of a. For a > 2, the system is 

ultradamped. For 2 > a >0, the system is hyperdamped, passive, and positive-real. For 0 > a > - -Jl , the 
system is underdamped and minimum phase. For a < - yfl , the system is unstable. 
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Phase, degrees 



Frequency, radians 

Figure 3.6.1. Bode frequency response for ultradamped systems. 



Real H(jw) 

Figure 3.6.2. Nyquist frequency response for ultradamped systems. 
The arrow shows the direction of increasing frequency. 
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normalized step reponse 



33 
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The underdamped active region bears further consideration. Specifically, the negative value of a 
indicates an active system, while the tv-plane poles remain in the underdamped region. The fact that there 
exist underdamped poles implies that this system is in resonance, and a resonant peak should appear in the 
system frequency response. Figure 3.7.3 shows the Bode frequency response for several values of a. 
A resonant peak is clear in this response. The phase response behavior is more interesting. For integer- 
order systems, the phase will decrease by 180 degrees at a resonance, independent of the amount of 
damping. This also occurs for fractional-order systems, however, as Figure 3.7.3 shows, the phase 
increases before the resonant frequency, then the phase decreases at its maximum rate at the resonant 
frequency, and finally increases again back to the high frequency asymptote, which is -90 degrees for this 
system. The phase decrease at a marginally stable resonance remains 180 degrees, although the phase 
decrease gets smaller with more damping in the system (see the a = - 1.4 and a = -1 cases in Figure 3.7.3). 
Thus, even though the high frequency asymptotes (as well as the Laplace transfer function) indicate that 
this system is only of first order, it can still go into resonance. Clearly, adding more fractional-order terms 
with a smaller value of q would allow even more resonances. Consequently, it appears that the highest 
power of the Laplace variable in a transfer function is no longer an indicator of the effective order of a 
fractional-order system, or of the number of resonances to expect in its frequency response. We will call 
resonances of the type considered here fractional resonances. 



Figure 3.7.3. Bode plot for the system given in Equation (3.7.1). 


3.8 Nyquist Plane Design 

Frequency response techniques have proven very useful for many control design problems. In 
control design, frequency response information is often presented in one of three equivalent forms; the 
Bode plot, the Nyquist plot, or the Nichols plot. Bode plots are widely used by engineers, and 
consequently, are the most widely understood presentation of the systems frequency response data. These 
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plot the magnitude and phase of the system frequency response against a logarithmically increasing 
frequency. The Nyquist plot is an equivalent representation, however, it is a plot of the real versus the 
imaginary part of the system frequency response. This plot contains much useful information for feedback 
systems, as is shown in the next section. The Nichols plot is a conformal mapping of the Nyquist plot, and 
is probably in less use due to the distortion of the information. The remainder of this section contains a 
brief review of the utility of the Nyquist plot, and its use in a fractional-order system setting. 


Yp(s) ¥k( s ) Vg (■*) 



Figure 3.8.1. Closed-loop control configuration. 


Referring to the single-input-single-output system of Figure (3.8.1), the signal zlx) is a desired 
reference input, d(s) is some generally unknown disturbance acting on the plant G(s) and can include the 
plant uncertainty, m(s) is the error incurred in measuring the plant output y(s), and u(s) is the control 
signal applied to the plant from the controller K(s ), and P(s) is a prefilter used to shape the system inputs. 
The functions i jf P (s), y/ K (s), and y/ G (5) are the initialization functions associated with the 
corresponding transfer functions in the block diagram. The closed-loop plant response can the be 
determined to be 


J(I) = 1 ?G^k lPMr(s,] + iT (s 'l 

(3.8.2) 


This equation is organized so that the effects of the time varying initializations can be directly understood. 
The prefilter initialization \ff P (s ) , is combined with the measurement noise term, and has the same 
effect as measurement noise. The plant initialization and the controller initialization are 
combined with the disturbance term, and have the same effect as plant disturbances or 
uncertainty. 
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Normally the following definitions are made; the loop gain 


L(s)=G(s)K(s ) 

is the product of the plant and controller transfer functions. The sensitivity function 

S(s) = ! 

1 + G(s)K(s) 


(3.8.3) 


(3.8.4) 


represents the sensitivity of the output to plant disturbances. The complementary sensitivity function 


G(s)K(s) 

1 + G(s)K(s) 


(3.8.5) 


is the sensitivity of the output to the reference input as well as measurement errors, and is normally 
known as the closed-loop transfer function. 

It is noted that 

7X$) + S(s) = 1 . (3.8.6) 

This equation indicates that there is a trade off between reducing the effects of plant disturbances and 
reducing the effects of measurement noise. The magnitude of T(s) is also an indication of the control 
effort required, u(t). This trade-off is usually managed by making T(jco) small at one set of frequencies 
and then making S(j(0 ) small at another set of frequencies. Most often, the measurement noise m(t) will 
have a large bandwidth, and thus T(ja>) will be made small at high frequencies. Correspondingly, the 
disturbances usually have a low frequency nature, and thus S(jco ) is made small at low frequencies. The 
trade-off is usually accomplished by forcing the loop gain, L(s) = G(s)K(s) , to be large at low 
frequencies and small at high frequencies. A consequence of not allowing T( j( 0 ) to be equal to one at all 
frequencies is that the system cannot track an input perfectly, and equivalently, the control actuators will 
not be required to be effective at unrealistically high frequencies. 

One of the most useful benefits of using frequency response techniques in linear control system 
design is that closed-loop information can be directly obtained from open-loop information. Referring to 
the Nyquist approach, not only can stability be determined by looking at encirclements of the minus-one 
point, but also the frequency response plots of T(s) and S(s) can be determined directly from the 
frequency response plot of L(s). To see these, lines are added to the Nyquist plane representing contours 
of constant magnitude of both T(s) and S(s). A detailed discussion of these contours can be found in 
D’Azzo and Houpis (1981) or Maciejowski (1989), and a brief discussion follows. For T(ja>), its 
magnitude is zero at the origin of the Nyquist plot, the magnitude is infinity at the minus-one point of the 
Nyquist plane, and the unity magnitude contour is the vertical line at Real( L(jco) )= -0.5. For S(j(o) , its 
magnitude is zero at a radius of infinity on the Nyquist plane, its magnitude is infinity at the minus-one 
point, and the unity magnitude contour is a circle of radius one centered at the minus-one point. The 
consequences of this are that the sensitivity to disturbances at low frequencies can be made zero by the 
addition of an integrator to the controller. The addition of the integrator will also allow perfect tracking of 
inputs at low frequencies as well. If the pole excess is greater than one, the Nyquist plot of L(s) will enter 
the unity sensitivity circle at high frequencies, indicating that the cost of lowering the sensitivity to low 
frequency disturbances is an increase in sensitivity to high frequency disturbances. 
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All of the frequency response tools from integer-order system theory carry over to fractional- 
order systems directly, as long as the frequency response plots are taken to be on the primary Riemann 
sheet. In the next section, some examples are presented showing a comparison of Nyquist plane and 
u’-plane root-locus techniques for several control configurations. 


3.9 Controllers in a Fractional-order Universe 

Let us consider a fractional-order plant with components that are multiples of 1 /2-order. To 
control such a system, we have at our disposal all the standard compensators from integer-order control 
systems, with the addition of 1/2-order terms such as semi-integrals, semi-derivatives, semi-leads, semi- 
lags, semi-pole placement, etc. This extra design freedom allows improved control designs. Usually an 
improved control design will increase the stability, the tracking performance, the speed of response, the 
insensitivity to disturbances, or some similar measure of goodness. 

Many generalizations of integer control design are possible with the freedom allowed by 
fractional-order systems. Some of the many generalizations are listed below. 

a) Integral control: 

Fractional integrals are now available for compensators. The interesting feature of fractional integrals is 
that they still allow closed-loop tracking of step reference signals, while allowing the freedom to tune the 
low and high frequency behavior by tuning the value of q , although the tracking will be slower 

b) Derivative Control: 

Although pure derivative control is seldom used, derivatives of any fractional-order are now available and 
these will have less noise amplification at high frequencies than integer-order derivatives. 

c) PI, PD, and PID Control: 

Fractional elements allow the use of any values of q for the integral and the derivative in these controllers. 
If fractional PED control is implemented, the fraction in the derivative need not be the same as the fraction 
in the integral. Podlubny (1999a and 1999b) discusses controllers of the form P/ p D v more thoroughly. 

d) Leads and Lags: 

Lead compensators are often used to help stabilize marginally unstable systems. Lag compensators are 
often used to help reduce the magnitude of the high frequency loop gain of the system. Using fractional- 
order components, it is now possible to design fractional leads and fractional lags. The benefit to these is 
that it is easier to shape the open-loop and closed-loop frequency responses using them than using 
exclusively integer-order elements, due to the extra freedom offered by the continuum of values of q. 

Another design criterion unique to fractional-order systems deals with time domain singularities 
occurring at time zero. If a transfer function, F(s ) , does not contain a term in the denominator with an 

exponent of at least 1 , that is s 1 , then by the initial value theorem 

/( 0) = lim sF(s) , (3.9.1) 


the impulse response, f(t) = L 1 {F(s)}, will have a singularity at time zero. This may or may not be 
desirable, however, using the appropriate fractional feedback can eliminate the singularity. 


NASA/TM— 2002-21 1377/REV 1 


21 


It can occur by design, or otherwise, that a fractional-order system with a given base order, say 
1/2, must be combined with a fractional-order system that has a different base order. That different order 
can be rationally related, say 1 or 1/3, or it can be non-rationally related, say \/n . In the first situation, all 
fractional components are related to one another by rational fractions. The plant may contain only integer- 
order devices, whereas the controller could be of half-order. The second situation deals with systems 
whose components are not rationally related fractions. It should be noted that although any of the methods 
discussed above can very closely approximate systems with irrationally related components by choosing a 
nearby rational fraction, they cannot solve it exactly as is discussed below. An example will illustrate the 
situation. Consider the plant 


G(s) = 

s + 1 


(3.9.2) 


to be compensated by the specific fractional-order integral 


H(s) = 



1 

0 . 318309886 ... ’ 

S 


(3.9.3) 


In this situation, one would normally choose the base fraction for constructing the appropriate w-plane. 
This value of q is usually the greatest common denominator of all the fractions appearing in the system. 
Unfortunately, for the above system in closed-loop, the greatest common denominator goes to zero, or 
some very small number, because pi is irrational. This would imply that the corresponding vv-plane would 
be almost entirely a hyperdamped region, with a small sliver representing the usual s-plane. This would 
imply that the corresponding branch points in the 5-plane are of infinite order, which means that the 
number of Riemann sheets below the primary sheet goes to infinity. Finally, a fractional vector space for 
such a system would require an infinite number of fractional vectors. This is not a very practical situation. 

The most straightforward approach to analyzing such systems is to use Nyquist design from the 
beginning. The fact that there are non-rationally related terms in the system is of no consequence when 
using Nyquist analysis, as the frequency response uses the response on the primary Riemann sheet. Thus 
the frequency response and the Nyquist analysis tools can be used as usual, without any confusion or 
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Figure 3.9.1. Frequency response analysis of the combination of the system of Equation (3.9.3) 
used to compensate the system of Equation (3.9.2). 
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modification. This is shown for the given system in Figure 3.9.1, with the inclusion of the unity 
sensitivity circle (center at s = -1, radius 1) and the unity complementary sensitivity line (at s = -0.5) on 
the Nyquist plot. Consideration of the Nyquist plot of the plant-controller combination with respect to the 
unity sensitivity circle and complementary sensitivity line shows that the closed-loop system will track 
steps, however sluggishly, and that there will be just a little bit of overshoot for large enough gains. 


3.10 Generalized Pi-Control and PID-Control 

This section first presents a generalization of the proportional-plus-integral-(PI)-controller. This 
is followed by a generalization of proportional-plus-integral-plus-derivative-(PID)-control. These 
generalizations are possible only because of the availability of fractional-order elements. The standard 
integer-order Pi-controller has a transfer function of the form 

H(s) = k p +^. (3.10.1) 

s 

This can be written 

H(s) = k p + k;S~ l . (3.10.2) 

This controller can now be generalized using fractional-order integrals. Assuming a base fraction q=l/N, 
for now, a generalized Pi-controller can be written 

H(s) = k 0 + + k 2 s- 2q + ••• + * w _ 1 s" w “ lw + k N s~ l (3.10.3) 


or 


H(s) = - N « = l - 

n=0 

Again performing some algebra, this can be written as 

hw - 

s 


(3.10.4) 


(3.10.5) 


or 


'tty’ 1 -'" 

H(s) = — , Nq = 1. (3.10.6) 

Clearly, this controller allows a much greater variety of possible compensation results. Inserting this into 
the standard closed-loop control configuration with plant G(s) , gives the closed-loop transfer function, 
which is C(s) = H(s)G(s )/{ 1 + H(s)G(s)), to be 
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" /V 

£ k„s (N - n) ‘ l G(s) 

C(s) = =L , Nq= 1. (3.10.7) 

s + G(s) 

_ o=0 

With C(j) = A^(5)/D(s)this becomes 

W(s) 

C(s) = i="~ u ^ = , A^ = l. (3.10.8) 

jD(j)+ X /t « 5< ' V ’'' ,? yV(j) 

.' 1=0 


This major generalization of Pi-control possesses considerably more design capability and freedom, as 
both closed-loop poles and closed-loop zeros can be placed by proper selection of the gains. 

The compensator can be further generalized by considering the powers of q to be unrelated. In 
this case the compensator is 

H(s) = k 0 + + k 2 s~ q 1 + ■■■ + k + k N s~ l , (3.10.9) 

with q N assumed to be unity and q 0 assumed to be zero, and 0 <q t <1, 0<i<N . This can also be 
rewritten as 

H(s) = — . (3.10.10) 

5 

This result effectively generalizes the work of Maia (1998) to the design of system controllers. 

The above approach can also be applied to PID-controllers, where fractional-order derivatives are 
allowed. A similar approach has been presented by Podlubny (1999b), and like most discussions of PID- 
control we will not consider the causality problem. The Pi-controller of Equations (3.10.3) and (3.10.4) 
can be generalized to include fractional derivatives as follows, 

H(s) = k N s +l + + --- + k 2 s +2q + k,s +q + k 0 + k,s~ q + k 2 s~ 2q + ••• + k^s^^ 9 + k N s~ l 

(3.10.11) 

or equivalently, with q = 1 / N , 

N 

H(s) = , Nq= 1. (3.10.12) 

n=-N 


NASA/TM — 2002-2 1 1 377/REV 1 


24 



3.11 Generalized Dynamic Compensator 


Based on the results of the previous section, it should be clear that a compensator, such as 
Equation (3.6.10), can be further generalized by expanding not only the numerator into a series of 
fractional operators, but also the denominator. Hence, a generalized dynamic compensator has the form 


His) = 




ci„s 


n IN 


n=0 


(3.11.1) 


if the powers of s are related and are bounded by one. Alternatively, allowing arbitrary powers of s yields. 


His) 




(3.11.2) 


where here the numerator and denominator can have different powers of s. In a closed-loop feedback 
configuration with plant Gis) = Nis) / Dis ) , the closed-loop transfer function would be 


C(s) = 



,n=0 

His) 

ia„s- 

_n= 0 

Dis) + 

' N 

,H=0 

Nis) 


(3.11.3) 


To determine the effectiveness of such a controller, the analysis must be performed in the Nyquist plane. 
Any desired Nyquist locus may be approximated by using a generalized compensator. The quality of the 
approximation is limited only by the acceptable number of terms. 


3.12 Continuum Feedback: Order-Distributions 

Proceeding from the generalized compensator of the last section, the idea of taking the sum to the 
limit occurs. In this case, the summations would be replaced by integrals over the variable q, which is the 
power of 5. A continuum-of-g compensator generalizing Equation (3.1 1.2) would then have the form 

b 

| K N iq)s q dq 

His) = (3.12.1) 

u 

J K D iq)s q dq 


where the functions K(q) must be chosen so that His) remains causal. Integrals of the form shown in 
Equation (3.12.1) are called “order-distributions” in Hartley and Lorenzo ( 1999). 
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The generalized transfer functions of the last section can be written as infinite series. Thus both a 
fractional power series 

Iks'- . 

His) = sf = 

n=0 

and a fractional asymptotic series 

S 

- 

= ^ = ^c n s-"" 

” =0 

11=0 

representation can be used. Clearly these concepts may be generalized to the continuum. Thus a power 
integral representation for the continuum compensator would be 

a 

H{s) = jK PS (q)s q dq, (3.12.4) 

0 

while the asymptotic integral representation would be 

H(s) = ^K AS iq)s- q dq\ (3.12.5) 

0 

or in a combined form 

a 

H(s) = jK(q)s q dq, (3.12.6) 


(3.12.2) 


(3.12.3) 


where the K(q) functions are chosen so that the integrals are convergent. For example, an order- 
distribution PID-controller generalizing Equation (3. 10.12) is 

i 

H(s)=jK(q)s q dq . (3.12.7) 

-i 

Clearly, given the K(q) functions, it is an easy task to perform analysis of a closed-loop system in the 
Nyquist plane. However, design of the Kiq) functions for a particular behavior is still an open question. 
Exact physical realization is also problematic, although accurate approximations are presently possible. 
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As another example of a possible situation, consider the plant to be controlled to be 


G(s) = -1— (3.12.8) 

s' +1 

which represents an undamped oscillator, or an undamped spring-mass system. A possible compensator 
using order-distributions would be 


H(s) = ^K(q)s q dq . (3.12.9) 

o 

If this compensator were placed before the plant, the resulting closed-loop system would have the transfer 
function 


| K(q)s q dq 

C(s) = ^ . (3.12.10) 

s' + J K(q)s q dq + \ 
o 

This compensator allows an infinite number of frequencies in the closed-loop system, and thus allows 
considerable freedom to design the appropriate K{q ) . The question arises as to how one chooses K(q). A 
straightforward approach is to minimize some desired error for a given input. The error transfer function 
for a step input is 


E(s) = - 

S 


s° + K(q)s q dq 


s 2 + K(q)s q dq +1 
One standard performance measure is the integral-of-the-squared-error (ISE) 
j ise = \E(j(0)\' d(0 = £°|e(/)|" dt 


(3.12.11) 


(3.12.12) 


Analytical solutions for order-distribution compensators are difficult to obtain. However, the ISE optimal 
spring-mass system containing all the half order terms has been approximated by numerical optimization 
to be 


g(s) = 


1 

1.4212s 1 ' 5 + 1.8419s - 0.1 113s 05 + 1 ’ 


(3.12.13) 


with J ISE = 0.7273. This performance measure is considerably smaller than J ISE = 1 , which is the 
known value obtained using a single integer order damper . Clearly, the effect of adding more fractional 
terms would allow a further improvement in performance. 
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4. Example 


4.1 System Definition 

In this section, a simple example is presented to demonstrate how to use the theory presented in 
this paper to obtain the solution to a physical fractional-order system. The system considered is the 
viscoelastically damped oscillator shown in Figure 4.1.1. 



Figure 4.1.1. Spring-mass-viscoelastic system. 

The input to the system is the force applied at the left, fit), and the output will be chosen to be the 
displacement of the mass, x(t ) . The force due to the spring is 

f s =k l x(t) (4.1.1) 

and the force due to the viscoelastic element is 

f v =k 2 0 D?x(r). (4.1.2) 

Summing forces yields 

m 0 D;x(t) = fit) - k [ x{t) - k 2 0 Dfx(t) , 
or 

m 0 Dfx{t) + k 2 0 D?x(t) + fc,.v(r) = fit ) . 

In terms of the initializations, this can be written 

m( 0 d,-x(t)+y/ix m ,2,a,0,t))+ k 2 ( o d?xit)+\//ix y ,q,a,0,t )) 

+ fc, ( 0 d° x(t)+y/ix s ,0,a,0,r)) 

where we have used -v m (r), x v (r), and x s it) as the individual displacements of the mass, viscoelastic, and 
spring elements during the time a<t< 0, while xit) is the displacement of the combined system after 


= fit), t> 0. 


(4.1.5) 


(4.1.3) 

(4.1.4) 
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time t = 0 . This is the fractional-order differential equation describing the behavior of the spring-mass- 
visco system, with input fit) for t > 0 , with proper initialization of the components. 

The general case of separate component histories is chosen to illustrate the most general 
initialization procedure for the solution of fractional differential equations. For this sample initialization, 
each component of the system has a separate history and will be treated separately, however it is required 
that all components have the same displacement at time zero and afterward. The spring is depressed to 
one unit displacement at time t = t s where a < t s < 0 , and held there until time zero using constant force 
f s . The mass is also displaced one unit at time / = t m where a < t m < 0 , and held there until time zero 
using force f „, , so that its velocity is zero. The viscoelastic element is displaced one unit at t = t v = a and 
held there until time zero so that ,v v (r)=l for a < t < 0 , using the force /,,(/) = k 2 (t-a )~ q /T(l -q) for 

a < t < 0 , which yields the initialization to be discussed below. It is understood that at t = 0 and 
afterward, all of the system elements are connected together, following which time, the various 
initialization effects superpose together to give the overall system initialization response, and only one 
external force, fit) , acts to give the forced response. 


The system equation (4.1.5), can be Laplace transformed as 

(/ns"x(5) + /7zL{vr(v m ,2,a,0,r)}) + (k 2 s q x(s) + k 2 L^f(x v ,q,a,0,t)^j + (k ] xis) + k l L\y/(x s ,0, a,0, f)})= f(s ) . 

(4.1.6) 


Solving for *(s) , gives 


f(s) , -m\i/ m (s)-k 2 y/ v (s)-k i \i/ s (s) 
ms~+k 2 s q +k ] ms'+k 2 s q +k i 

with the transformed initialization functions defined to be 

V'mCs) = L{vz(x m ,2,a,0,/)}= l{- x m i0) 0 J'5(r)-L m (0)5(r)} 

= -«»(0)-i m f0) 

V s (s) = Lfo/(x s ,0,a,0,?)}= 0, 
and 

\f/ y {s) = L$f/(x v , q, a,0, f)} = l\^~ f ^ v,, (r ) dr 1 . 

I dt{Ti\-q) j 


(4.1.7) 


(4.1.8) 

(4.1.9) 

(4.1.10) 


The first term in Equation 4.1.7 is the forced response, and the second term is the overall system 
initialization response. 

The viscoelastic initialization function is now evaluated. Again, with x v (t) = 1 for a < t < 0 , 


y/ v (s) = L^f{\,q,a,0,t)} = 


kf iziT! 

\dt{ ra -q) 



(t — ay q -t~ q 

m-^) 


(4.1.11) 
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or 


¥M*) 


e~ as T(\- q-as) 
HI-*) 


(4.1.12) 


which can be found in Oberhettinger and Badii, page 21. Thus, inserting Equations (4.1.8), (4.1.9), and 
(4.1.12) into Equation (4.1.7), with .v(0) = l and .v(0) = 0, the overall response function in the Laplace 
domain is 


,v(s) = 


f(s) 
ms~ +k-,s q + k. 


ms + k-> s q 1 


1 - 


e °T(1 - q -as) 


T(\-q) 


ms 2 +k-,s q + £, 


(4.1.13) 


For demonstration purposes, let m = 1, k t = 1, q = 0.5, k 2 = 1 . Also let the forcing input be f (t) = l, t > 0 , 

e- a T(0.5 -as) X 


which has the Laplace transform f(s) = \/ . Then 


s + s 


-0.5 


-V(5) = 


1 - 


no.5) 


s(s~ +5° 5 + 1) 


- i 0.5 , i 
5+5 +1 


(4.1.14) 


The w-plane poles of this system are w=0.7271± j0.9341,-0.7271± j0.4300, all of which are stable as 
they lie to the left of the 45-degree instability wedge in the vv-plane. The first pair of poles are 
underdamped which would imply some oscillation in the response, while the second pair of poles are 
hyperdamped. 

4.2 Forced Response 

The forced response part of Equation (4.1.14) is found as the inverse Laplace transform of 


x f (5) = 


s(s~ +5°' 5 + 1) 


which is easily evaluated as the inverse w-transform of 

1 


x f (w) = 


6 , 3 . : 

w + w + vv 


(4.2.1) 


(4.2.2) 


where the substitution w = 5° 5 has been used in the first term on the right of Equation (4. 1 . 1 1 ). This is 
expanded into partial fractions as 


-1 J_ - 0.0291 + j0. 1210 -0,0291- j0. 12 10 

V vv 2 w -0.7271 - j0.9341 w -0.7271 + j0.9341 


0.5291- j0.0440 0.5291 + j0.0440 

w + 0.727 1 - j0.4300 w + 0.727 1 + j0.4300 ' 


(4.2.3) 
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Replacing vv by s 0 ' 5 , we can inverse Laplace transform term-by-term to give the forced time response 

, , 1 

xAt) ■=+ 1 + 

r(0.5)Vr 

+ (- 0.029 1 + jO. 1 210)F,/ [ 0.727 1 + j0.934 1 , t]+ (- 0.029 1 - jO. 1 210 )F X/ [ 0.727 1 - j0.934 1 , t] 

/: /z 

+ (0.529 l-j0.0440)F 1/ [ -0.7271 + j0.4300,r]+ (0.5291 + j0.0440)/v [-0.7271 - j0.4300,r], (4.2.4) 

/z 7z 

This is shown graphically in Figure 4.2.1, and the F-function was defined in Equation (2.1.21). It should 
be noted that this time response gives a steady state displacement of one for t large. 



Figure 4.2.1. Forced response of spring-mass-visco system. 


4.3 Initialization Response 

The initialization response part of Equation (4. 1. 14) is found as the inverse Laplace transform of 

e" a T(0.5 -as) ' 


s + s 


- 0.5 


= 


1 


no.5) 


s z +s 05 +l 


(4.3.1) 


For this discussion, we will let the initialization period begin at a = — «> for simplicity, which drives the 
term containing the gamma functions to zero (other initializations could just as well be considered). Then 


5 + 5 

*fc(5) = — 


- 0.5 


1.5 , , 
5 + 1 


5-+5 U5 +l 5- S +5'+5 l 


(4.3.2) 
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This equation can be expanded into vv-domain partial fractions. Then, replacing w with s 0 ' 5 , it can be 
inverse Laplace transformed term by term; 




w + 1 

S ^ \ 

w + w~ + vv 


1 0.0291- jO. 1210 0.0291 + j0.1210 

— + + 

vr w- 0.727 l-j0.9341 w - 0.7271 + j 0.9341 

- 0.529 l + j0.0440 - 0.5291 - j0.0440 

+ - + - , 

w + 0.727 1 - j 0.4300 w + 0.727 1 + j 0.4300 


(4.3.3) 


(4.3.4) 


.v,,(f) = 


T(0.5)Vr 


+ (0.029 1 - jO. 1 2 1 0 )F U [ 0.727 1 + j0.934 1 , /]+ (0.029 1 + jO. 1 2 1 0)F t/ [ 0.727 1 - j0.934 1 , t] 

Zi Zi 

+ (- 0.529 1 + j0.0440)F 1/ [ - 0.727 1 + j0.4300, /] + (- 0.529 1 - j0.0440)F,/ [ - 0.727 1 - j0.4300, r] . 


(4.3.5) 


This is shown graphically in Figure 4.3.1. 



Figure 4.3.1. Initialization response of spring-mass-visco system. 
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The total response is 

x(t) = x f (t) + x k (t), (4.3.6) 

which for this particular input and initialization is simply x(t) = 1 , t > 0 . 

4.4 Vector Space Solution 

This section presents the solution to the spring-mass-visco problem in the time domain using the 
vector space approach. To have a vector space of fractional dynamic variables, we will reduce all of the 
fractional derivative relationships to fractional derivatives based on the largest common order. In this 
case, the largest common order is q= 1/2. To define the fractional dynamic variable vector, we will 
represent the system of Equation (4. 1.5) in a block diagram representation. This is done so that it is easy 
to maintain the initializations with the correct vector variables. To accomplish this. Equation (4.1.5), with 
y/ s (t )=0 , is integrated twice to give 

.v(/) = 0 d; 2 |-^-(/(r) - k 2 ( Q d°- 5 x(t)+y/ v (r))- £, 0 ^,°v(/))- t// m (r)| . (4.4.1) 

This equation can be represented in a simple block diagram with x(t) as the output of a summation of 
those terms on the right side as shown in Figure 4.4.1. This integral form for representing differential 
equations has often been used for analog computing as it naturally incorporates the initialization effects of 
past inputs. It is unfortunate that this form is not often seen in the analytic solution of differential 
equations, which has consequently required the appearance of delta functions and their derivatives in the 
initialization functions (as in Equation (4.1.8)). Performing standard block diagram manipulations in 
Figure 4.4.1 can yield the block diagram of Figure 4.4.2, which is used for the structure of our vector 
space representation. 



Figure 4.4.1. Integral-form block diagram representation of the mass-spring-visco system. 
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Figure 4.4.2. A modified block diagram representation of the mass-spring-visco 
system useful for vector space realization.. 
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This vector space representation can now be Laplace transformed as in Equation (2.1.19) to give 


X(s) 


1.5 . t 0.5 1 

5+1 5 5 1 


'O' 


0 T 

1 1.5 0.5 

- 1 5 5 5 


0 


0 

n e / n ^ \ 1 c 

5 


U(s ) + 


- s - (s + 1 J s s 


0 


0 

-s -5 as (r a5 +i) -(s 05 + 1 ) s l5 _ 


1 


-yr v (s)-yr m (s) J 


(4.4.5) 


where, y/ v (s)=\(/(x v ,0.5,-«>,0,r) = — -=■ , and y/ m (s) = - sx(0) - x(0) which represents the initial 

V-s 

position and velocity effects of the mass on its position in the future. The Laplace transform of the forced 
vector response, X_ F (5) , is the first term of this equation, which contains the U(s) = term, 


X f (s) 


1 



+ 5 0 ' 5 + 1 


1.5 


U(s) , 


(4.4.6) 


and the Laplace transform of the vector initialization response, X_ ic (s), is the second term of this 
equation. 


K ic (s) = - 


S- + 5 05 + 1 


1 

0.5 


5 + 


V7 


(4.4.7) 


It can be seen that the first element in this 4-vector is identical to the displacement initialization response 
in Equation (4.3.1). Using the output equation (as in Equation (2.1.20)), the system output is given by the 
sum of the first (top) elements of the vectors in Equations (4.4.6) and (4.4.7), which is identical to the 
sum of Equations (4.2.1) and (4.3.2). These expressions can also be evaluated for any other specific 
inputs and initialization functions. 

It should be noted that the availability of the fractional vector formulation allows many choices for 
the basis <7 , and thus the total number of fractional dynamic variables can increase if the basis q is 

chosen to be smaller. In the example of this section, choosing <7 = 1/4, rather than 1/2, would yield eight 
fractional dynamic variables rather than three. It is also important to remember that the least number of 
fractional dynamic variables is obtained by choosing the basis q as the largest common fraction of the 

differential orders. For example, a system with two elements of order q l =1/2, and q 2 =1/3, would 
require a basis q of at most <7 = 1/6. 

The next section will use an input-output approach for controlling the position of the mass in this 
example, by using a force input and a position output. After that, a vector space approach is used for 
controlling the position of the mass using state feedback. 


NAS A/TM— 2002-2 1 1 377/REV 1 


35 


4.5 Input-output Control 

An input-output control approach is taken first, followed by the vector space approach. Using the 
input-output approach, there are two tools that can be used for analysis and design, the vv-plane root- 
locus, and the frequency response. The system transfer function is repeated below, where it should be 
noted that the initialization responses show up as plant disturbances, which always decay to zero for 
0<q< 1 , 


-v(s) = 


(s 2 +$ as + l) 


f(s) + 


s.y(O) + i(0) - t/A(s) 
s 2 +s 0 ' 5 +1 


(4.5.1) 


The open-loop and closed-loop system behavior is studied in Figure 4.5.1 for both unity feedback control 
and the specific semi-PID control given here 


f(s) = ~ 


5.5 (s - 5°' 5 + 0.8) 
s 05 (5 0 - 5 + 1) 


x(s) + r(s ) , 


(4.5.2) 


which was obtained by trial and error, and where r(s) is the Laplace transform of a reference signal. Note 
that the semi-integral in the compensator forces the closed-loop system to track step reference signals, as 
is the case for integer integral feedback. 


4.6 Vector Space Control 

For the vector space controller, there are also two basic approaches, vector feedback and output 
feedback. The vector feedback situation requires sensors on all of the system vector variables so that each 
of these variables can be fed back to the system input, thereby moving the system poles in the vv-plane. 
The output feedback controller requires a system observer that will recreate the system vector, which is 
then fed back to the system input to move the vv-plane poles. With both of these approaches, an integrator 
is often included to act on the difference between the reference input and the system output, thereby 
forcing the system to track steps in the reference input. 

State feedback is now considered. The system equations are 


) d]' 2 x(t) = 


“0 

1 

0 

o' 


0 


0 

0 

0 

1 

0 

x(t) + 

0 

u{t) + 

0 

0 

0 

0 

1 

0 

0 

-1 

-1 

0 

0 


1 




y{t) = [l 0 0 0] x(t) + [o] u(t) , 


(4.6.1) 


(4.6.2) 


where y/ r (t) represents the initialization function related to the viscoelastic term and t y m (() represents 
the initialization related to the inertial term. Typically, vector feedback is implemented in the form 

u(t) = - K x(t) + r(t) 
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Figure 4.5.1. Spring-mass-visco system using the semi-PID controller from the text. 


where r(t) is some reference input, and K is the feedback gain matrix to be determined. Using Equation 
(2.3.3), the closed-loop system becomes 

0 d? x(t) = [A - jBA'Jxfr) - y/(x,q,a,0,t) + B r(t) , y(t) = [C - DK]x(t) + D r(t) (4.6.3) 

By choosing K appropriately and using any of the standard pole placement tools, such as the Ackerman 
method, it is possible to place the system eigenvalues anywhere in the vv-plane. Several closed-loop 
responses are shown in Figure 4.6.1 with the initialization functions set to zero (ISE Optimal from 
Equation (3. 12. 13)). Notice that the ISE optimal has a fast rise time, but a sluggish settling time. 


5. Discussion 

An overview of control system design techniques for fractional-order systems have been 
presented in this paper. Some of the important features of the discussion are listed below. 

a) Fractional-order behavior occurs in viscoelastic systems, electrochemical systems, heat 
transfer problems, and in other areas. 

b) Control of these systems is very important and becomes more so when it is recognized that 
fractional-order systems are a generalization that includes the normal integer-order systems 
as a special case. 
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Figure 4.6.1. Time responses of several controllers for the spring-mass- visco system. 


c) The issue of control is complicated by two factors. One is the presence of the time- varying 
initialization term (the fading memory), and the other is the presence of non-integer-order 
derivatives and integrals. 

d) Vector space control design is done in almost the same way as for integer-order systems, 
with the major differences being that the system equations are fractional-order differential 
equations, and that the pole placement must be done in the w-plane corresponding to the 
basis value of q. 

e) Frequency response techniques are valid for the analysis of fractional-order systems, as long 
as the principal values of the fractional terms are used. 

f) With the validity of the frequency response, both Bode and Nyquist design techniques can be 
used with little modification. 

g) Root locus design methods can be applied directly to fractional-order systems where the root 
locus is performed in the w-plane. 

It is clear that the presence of fractional-order systems, and the availability of fractional-order 
compensators, allows a wide variety of new control approaches to be used. It will be interesting to see the 
effect of these fractional-order systems in future control applications. One should remember, however, 
that when using fractional-order components, it is important to properly include the effects of component 
initialization. 
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